{
 "cells": [
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T03:12:13.172843Z",
     "start_time": "2025-07-31T03:12:13.168735Z"
    }
   },
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy.optimize import least_squares\n",
    "from scipy.linalg import qr, svd, solve_triangular\n",
    "import warnings\n",
    "\n",
    "path = \"D:/7_30data/\""
   ],
   "id": "ae882b1b7ef98b7f",
   "outputs": [],
   "execution_count": 5
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T03:12:13.238639Z",
     "start_time": "2025-07-31T03:12:13.194090Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def high_precision_solve(V, L_rho, method='svd_truncated', **kwargs):\n",
    "    \"\"\"\n",
    "    高精度求解线性系统 S @ a = b，其中 S = V.conj().T @ V, b = V.conj().T @ L_rho\n",
    "\n",
    "    Parameters:\n",
    "    -----------\n",
    "    V : ndarray, shape (m, n)\n",
    "        输入矩阵\n",
    "    L_rho : ndarray, shape (m,)\n",
    "        右端向量\n",
    "    method : str\n",
    "        求解方法选择：\n",
    "        - 'svd_truncated': 截断SVD (推荐)\n",
    "        - 'qr_pivoting': QR分解with列主元\n",
    "        - 'iterative_refinement': 迭代精化\n",
    "        - 'tikhonov_adaptive': 自适应Tikhonov正则化\n",
    "        - 'multiple_precision': 多精度计算\n",
    "\n",
    "    Returns:\n",
    "    --------\n",
    "    a : ndarray\n",
    "        解向量\n",
    "    \"\"\"\n",
    "\n",
    "    # 计算正规方程矩阵\n",
    "    S = V.conj().T @ V\n",
    "    b = V.conj().T @ L_rho\n",
    "\n",
    "    # 输出基本信息\n",
    "    rank_S = np.linalg.matrix_rank(S)\n",
    "    augmented = np.hstack((S, b.reshape(-1, 1)))\n",
    "    rank_augmented = np.linalg.matrix_rank(augmented)\n",
    "    cond_S = np.linalg.cond(S)\n",
    "\n",
    "    print(f\"rank(S) = {rank_S}\")\n",
    "    print(f\"rank([S | b]) = {rank_augmented}\")\n",
    "    print(f\"Condition number of S: {cond_S:.4e}\")\n",
    "\n",
    "    if method == 'svd_truncated':\n",
    "        result = svd_truncated_solve(S, b, **kwargs)\n",
    "    elif method == 'qr_pivoting':\n",
    "        result = qr_pivoting_solve(V, L_rho, **kwargs)\n",
    "    elif method == 'iterative_refinement':\n",
    "        result = iterative_refinement_solve(S, b, **kwargs)\n",
    "    elif method == 'tikhonov_adaptive':\n",
    "        result = tikhonov_adaptive_solve(S, b, **kwargs)\n",
    "    elif method == 'multiple_precision':\n",
    "        result = multiple_precision_solve(S, b, **kwargs)\n",
    "    else:\n",
    "        raise ValueError(f\"Unknown method: {method}\")\n",
    "\n",
    "    # 根据问题类型决定是否取实部\n",
    "    # 如果原始问题期望实数解，则取实部；否则保持复数\n",
    "    if np.allclose(np.imag(result), 0, atol=1e-10):\n",
    "        return np.real(result)\n",
    "    else:\n",
    "        print(\"Warning: Solution has significant imaginary components\")\n",
    "        print(f\"Max imaginary part: {np.max(np.abs(np.imag(result))):.4e}\")\n",
    "        return result\n",
    "\n",
    "def svd_truncated_solve(S, b, tolerance=1e-12, verbose=True):\n",
    "    \"\"\"\n",
    "    使用截断SVD求解，这是处理病态矩阵最稳定的方法\n",
    "    修复复数矩阵处理\n",
    "    \"\"\"\n",
    "    print(\"\\n=== SVD Truncated Solve ===\")\n",
    "\n",
    "    # SVD分解\n",
    "    U, s, Vh = svd(S, full_matrices=False)\n",
    "\n",
    "    # 找到有效的奇异值\n",
    "    s_max = np.max(s)\n",
    "    valid_indices = s > tolerance * s_max\n",
    "    n_valid = np.sum(valid_indices)\n",
    "\n",
    "    print(f\"Total singular values: {len(s)}\")\n",
    "    print(f\"Valid singular values (> {tolerance:.2e} * σ_max): {n_valid}\")\n",
    "    if n_valid > 0:\n",
    "        print(f\"Smallest valid singular value: {s[valid_indices][-1]:.4e}\")\n",
    "    print(f\"Largest singular value: {s_max:.4e}\")\n",
    "    print(f\"Condition number (before truncation): {s_max/s[-1]:.4e}\")\n",
    "    if n_valid > 1:\n",
    "        print(f\"Condition number (after truncation): {s_max/s[valid_indices][-1]:.4e}\")\n",
    "\n",
    "    # 截断SVD求解 - 正确处理复数\n",
    "    s_inv = np.zeros_like(s)\n",
    "    s_inv[valid_indices] = 1.0 / s[valid_indices]\n",
    "\n",
    "    # a = V @ diag(1/s) @ U^H @ b (注意共轭转置)\n",
    "    a = Vh.T.conj() @ (s_inv * (U.T.conj() @ b))\n",
    "\n",
    "    # 计算残差\n",
    "    residual = np.linalg.norm(S @ a - b)\n",
    "    print(f\"Residual norm for S @ a - b: {residual:.4e}\")\n",
    "\n",
    "    return a  # 保持复数形式\n",
    "\n",
    "def qr_pivoting_solve(V, L_rho, tolerance=1e-12, verbose=True):\n",
    "    \"\"\"\n",
    "    使用QR分解with列主元直接求解 V @ a = L_rho\n",
    "    避免计算 S = V^H @ V，减少条件数恶化\n",
    "    修复复数矩阵处理\n",
    "    \"\"\"\n",
    "    print(\"\\n=== QR Pivoting Solve ===\")\n",
    "\n",
    "    # QR分解with列主元 - 对V^H进行QR分解\n",
    "    Q, R, P = qr(V.T.conj(), pivoting=True)  # 注意共轭转置\n",
    "\n",
    "    # 找到有效的对角元素\n",
    "    diag_R = np.abs(np.diag(R))\n",
    "    max_diag = np.max(diag_R)\n",
    "    valid_indices = diag_R > tolerance * max_diag\n",
    "    n_valid = np.sum(valid_indices)\n",
    "\n",
    "    print(f\"Matrix rank from QR: {n_valid}\")\n",
    "    if n_valid > 0:\n",
    "        print(f\"Smallest valid diagonal element: {diag_R[valid_indices][-1]:.4e}\")\n",
    "    print(f\"Largest diagonal element: {max_diag:.4e}\")\n",
    "\n",
    "    if n_valid == 0:\n",
    "        print(\"Warning: No valid diagonal elements found!\")\n",
    "        return np.zeros(V.shape[1], dtype=V.dtype)\n",
    "\n",
    "    # 求解 R[:n_valid, :n_valid] @ y = Q[:, :n_valid]^H @ L_rho\n",
    "    Q_valid = Q[:, :n_valid]\n",
    "    R_valid = R[:n_valid, :n_valid]\n",
    "\n",
    "    y = solve_triangular(R_valid, Q_valid.T.conj() @ L_rho)  # 注意共轭转置\n",
    "\n",
    "    # 构造完整解\n",
    "    a_permuted = np.zeros(V.shape[1], dtype=V.dtype)\n",
    "    a_permuted[:n_valid] = y\n",
    "\n",
    "    # 恢复原始顺序\n",
    "    a = np.zeros_like(a_permuted)\n",
    "    a[P] = a_permuted\n",
    "\n",
    "    # 计算残差\n",
    "    residual = np.linalg.norm(V @ a - L_rho)\n",
    "    print(f\"Residual norm for V @ a - L_rho: {residual:.4e}\")\n",
    "\n",
    "    # 同时计算正规方程的残差用于比较\n",
    "    S = V.conj().T @ V\n",
    "    b = V.conj().T @ L_rho\n",
    "    normal_residual = np.linalg.norm(S @ a - b)\n",
    "    print(f\"Residual norm for S @ a - b: {normal_residual:.4e}\")\n",
    "\n",
    "    return a  # 保持复数形式\n",
    "\n",
    "def iterative_refinement_solve(S, b, max_iter=10, tolerance=1e-14, verbose=True):\n",
    "    \"\"\"\n",
    "    迭代精化求解法\n",
    "    \"\"\"\n",
    "    print(\"\\n=== Iterative Refinement Solve ===\")\n",
    "\n",
    "    # 初始解（使用最小二乘）\n",
    "    a = sp.linalg.lstsq(S, b, lapack_driver='gelsy')[0]\n",
    "\n",
    "    print(\"Initial solution computed\")\n",
    "\n",
    "    for i in range(max_iter):\n",
    "        # 计算残差\n",
    "        r = b - S @ a\n",
    "        residual_norm = np.linalg.norm(r)\n",
    "\n",
    "        print(f\"Iteration {i+1}: Residual norm = {residual_norm:.4e}\")\n",
    "\n",
    "        if residual_norm < tolerance:\n",
    "            break\n",
    "\n",
    "        # 求解修正方程 S @ da = r\n",
    "        try:\n",
    "            da = sp.linalg.lstsq(S, r, lapack_driver='gelsy')[0]\n",
    "            a = a + da\n",
    "        except:\n",
    "            print(\"Iterative refinement failed, using current solution\")\n",
    "            break\n",
    "\n",
    "    final_residual = np.linalg.norm(S @ a - b)\n",
    "    print(f\"Final residual norm for S @ a - b: {final_residual:.4e}\")\n",
    "\n",
    "    return np.real(a)\n",
    "\n",
    "def tikhonov_adaptive_solve(S, b, verbose=True):\n",
    "    \"\"\"\n",
    "    自适应Tikhonov正则化 - 修复复数处理\n",
    "    \"\"\"\n",
    "    print(\"\\n=== Adaptive Tikhonov Regularization ===\")\n",
    "\n",
    "    # 使用L-curve方法选择正则化参数\n",
    "    U, s, Vh = svd(S, full_matrices=False)\n",
    "\n",
    "    # 候选正则化参数\n",
    "    lambda_values = np.logspace(-15, -1, 50)\n",
    "    residual_norms = []\n",
    "    solution_norms = []\n",
    "\n",
    "    for lam in lambda_values:\n",
    "        # 正则化求解 - 修正复数处理\n",
    "        s_reg = s / (s**2 + lam**2)\n",
    "        a_reg = Vh.T.conj() @ (s_reg * (U.T.conj() @ b))\n",
    "\n",
    "        residual_norms.append(np.linalg.norm(S @ a_reg - b))\n",
    "        solution_norms.append(np.linalg.norm(a_reg))\n",
    "\n",
    "    # 选择L-curve拐点处的参数 - 增加数值稳定性检查\n",
    "    curvatures = []\n",
    "    for i in range(1, len(lambda_values)-1):\n",
    "        # 添加数值稳定性检查\n",
    "        if (residual_norms[i] <= 0 or solution_norms[i] <= 0 or\n",
    "            residual_norms[i-1] <= 0 or solution_norms[i-1] <= 0 or\n",
    "            residual_norms[i+1] <= 0 or solution_norms[i+1] <= 0):\n",
    "            curvatures.append(0)\n",
    "            continue\n",
    "\n",
    "        # 计算曲率（简化版本）\n",
    "        try:\n",
    "            dx1 = np.log(solution_norms[i]) - np.log(solution_norms[i-1])\n",
    "            dy1 = np.log(residual_norms[i]) - np.log(residual_norms[i-1])\n",
    "            dx2 = np.log(solution_norms[i+1]) - np.log(solution_norms[i])\n",
    "            dy2 = np.log(residual_norms[i+1]) - np.log(residual_norms[i])\n",
    "\n",
    "            denominator = (dx1**2 + dy1**2)**1.5\n",
    "            if denominator < 1e-15:  # 避免除零\n",
    "                curvatures.append(0)\n",
    "            else:\n",
    "                curvature = abs(dx1*dy2 - dx2*dy1) / denominator\n",
    "                curvatures.append(curvature)\n",
    "        except (ValueError, RuntimeWarning):\n",
    "            curvatures.append(0)\n",
    "\n",
    "    if len(curvatures) == 0 or max(curvatures) == 0:\n",
    "        # 如果L-curve失败，使用GCV准则\n",
    "        print(\"L-curve method failed, using GCV criterion\")\n",
    "        return tikhonov_gcv_solve(S, b, U, s, Vh)\n",
    "\n",
    "    optimal_idx = np.argmax(curvatures) + 1\n",
    "    optimal_lambda = lambda_values[optimal_idx]\n",
    "\n",
    "    print(f\"Optimal regularization parameter: {optimal_lambda:.4e}\")\n",
    "\n",
    "    # 使用最优参数求解\n",
    "    s_optimal = s / (s**2 + optimal_lambda**2)\n",
    "    a_optimal = Vh.T.conj() @ (s_optimal * (U.T.conj() @ b))\n",
    "\n",
    "    residual = np.linalg.norm(S @ a_optimal - b)\n",
    "    print(f\"Residual norm for S @ a - b: {residual:.4e}\")\n",
    "\n",
    "    return a_optimal  # 保持复数结果，让调用者决定是否取实部\n",
    "\n",
    "def multiple_precision_solve(S, b, verbose=True):\n",
    "    \"\"\"\n",
    "    多精度计算（使用更高精度的数据类型）\n",
    "    \"\"\"\n",
    "    print(\"\\n=== Multiple Precision Solve ===\")\n",
    "\n",
    "    # 转换为更高精度\n",
    "    S_hp = S.astype(np.complex256) if hasattr(np, 'complex256') else S.astype(np.complex128)\n",
    "    b_hp = b.astype(np.complex256) if hasattr(np, 'complex256') else b.astype(np.complex128)\n",
    "\n",
    "    print(f\"Using precision: {S_hp.dtype}\")\n",
    "\n",
    "    # SVD求解\n",
    "    try:\n",
    "        U, s, Vh = svd(S_hp, full_matrices=False)\n",
    "\n",
    "        # 截断小的奇异值\n",
    "        tolerance = 1e-15\n",
    "        s_max = np.max(s)\n",
    "        valid_indices = s > tolerance * s_max\n",
    "\n",
    "        s_inv = np.zeros_like(s)\n",
    "        s_inv[valid_indices] = 1.0 / s[valid_indices]\n",
    "\n",
    "        a_hp = Vh.T @ (s_inv * (U.T @ b_hp))\n",
    "\n",
    "        # 转换回标准精度\n",
    "        a = a_hp.astype(np.complex128)\n",
    "\n",
    "        residual = np.linalg.norm(S @ a - b)\n",
    "        print(f\"Residual norm for S @ a - b: {residual:.4e}\")\n",
    "\n",
    "        return np.real(a)\n",
    "\n",
    "    except Exception as e:\n",
    "        print(f\"High precision solve failed: {e}\")\n",
    "        return svd_truncated_solve(S, b, tolerance=1e-12)\n",
    "\n",
    "def levenberg_marquardt_solve(V, L_rho, lambda_init=1e-3, max_iter=100, tolerance=1e-12):\n",
    "    \"\"\"\n",
    "    Levenberg-Marquardt方法求解非线性最小二乘\n",
    "    修复复数处理 - 将复数问题转换为实数问题\n",
    "    \"\"\"\n",
    "    print(\"\\n=== Levenberg-Marquardt Solve ===\")\n",
    "\n",
    "    # 将复数问题转换为实数问题\n",
    "    # [Re(V)  -Im(V)] [Re(a)]   [Re(L_rho)]\n",
    "    # [Im(V)   Re(V)] [Im(a)] = [Im(L_rho)]\n",
    "\n",
    "    V_real = np.vstack([np.real(V), np.imag(V)])\n",
    "    V_imag = np.vstack([-np.imag(V), np.real(V)])\n",
    "    V_expanded = np.hstack([V_real, V_imag])\n",
    "\n",
    "    L_rho_expanded = np.concatenate([np.real(L_rho), np.imag(L_rho)])\n",
    "\n",
    "    def residual_func(a_expanded):\n",
    "        return V_expanded @ a_expanded - L_rho_expanded\n",
    "\n",
    "    # 初始猜测\n",
    "    a0 = np.zeros(V_expanded.shape[1], dtype=np.float64)\n",
    "\n",
    "    try:\n",
    "        with warnings.catch_warnings():\n",
    "            warnings.simplefilter(\"ignore\", ComplexWarning)\n",
    "            result = least_squares(residual_func, a0, method='lm',\n",
    "                                 ftol=tolerance, xtol=tolerance, max_nfev=max_iter)\n",
    "\n",
    "        # 转换回复数形式\n",
    "        n_vars = V.shape[1]\n",
    "        a_real = result.x[:n_vars]\n",
    "        a_imag = result.x[n_vars:]\n",
    "        a_complex = a_real + 1j * a_imag\n",
    "\n",
    "        residual = np.linalg.norm(V @ a_complex - L_rho)\n",
    "        print(f\"Residual norm for Levenberg-Marquardt: {residual:.4e}\")\n",
    "\n",
    "        return a_complex\n",
    "\n",
    "    except Exception as e:\n",
    "        print(f\"Levenberg-Marquardt failed: {e}\")\n",
    "        return np.zeros(V.shape[1], dtype=complex)\n",
    "\n",
    "# 添加GCV方法作为备选\n",
    "def tikhonov_gcv_solve(S, b, U=None, s=None, Vh=None):\n",
    "    \"\"\"\n",
    "    使用GCV (Generalized Cross Validation) 准则选择正则化参数\n",
    "    \"\"\"\n",
    "    if U is None or s is None or Vh is None:\n",
    "        U, s, Vh = svd(S, full_matrices=False)\n",
    "\n",
    "    def gcv_function(alpha):\n",
    "        # 计算GCV函数值\n",
    "        filter_factors = s**2 / (s**2 + alpha**2)\n",
    "        numerator = np.sum((1 - filter_factors)**2 * (U.T.conj() @ b)**2)\n",
    "        denominator = (np.sum(1 - filter_factors))**2\n",
    "        return numerator / denominator\n",
    "\n",
    "    # 在对数空间搜索最优参数\n",
    "    alpha_values = np.logspace(-15, -1, 100)\n",
    "    gcv_values = [gcv_function(alpha) for alpha in alpha_values]\n",
    "\n",
    "    optimal_alpha = alpha_values[np.argmin(gcv_values)]\n",
    "    print(f\"Optimal regularization parameter (GCV): {optimal_alpha:.4e}\")\n",
    "\n",
    "    # 使用最优参数求解\n",
    "    s_optimal = s**2 / (s**2 + optimal_alpha**2)\n",
    "    a_optimal = Vh.T.conj() @ (s_optimal * (U.T.conj() @ b))\n",
    "\n",
    "    residual = np.linalg.norm(S @ a_optimal - b)\n",
    "    print(f\"Residual norm for S @ a - b (GCV): {residual:.4e}\")\n",
    "\n",
    "    return a_optimal"
   ],
   "id": "5c56296c07fa8747",
   "outputs": [],
   "execution_count": 6
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T03:12:13.260245Z",
     "start_time": "2025-07-31T03:12:13.254818Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# 示例使用函数\n",
    "def compare_all_methods(V, L_rho):\n",
    "    \"\"\"\n",
    "    比较所有方法的性能\n",
    "    \"\"\"\n",
    "    print(\"=\"*60)\n",
    "    print(\"COMPARING ALL HIGH-PRECISION METHODS\")\n",
    "    print(\"=\"*60)\n",
    "\n",
    "    methods = [\n",
    "        ('svd_truncated', {}),\n",
    "        ('qr_pivoting', {}),\n",
    "        ('iterative_refinement', {}),\n",
    "        ('tikhonov_adaptive', {}),\n",
    "        ('multiple_precision', {})\n",
    "    ]\n",
    "\n",
    "    results = {}\n",
    "\n",
    "    for method_name, kwargs in methods:\n",
    "        try:\n",
    "            print(f\"\\n{'='*20} {method_name.upper()} {'='*20}\")\n",
    "            results[method_name] = high_precision_solve(V, L_rho, method=method_name, **kwargs)\n",
    "        except Exception as e:\n",
    "            print(f\"Method {method_name} failed: {e}\")\n",
    "            results[method_name] = None\n",
    "\n",
    "    # 也测试 Levenberg-Marquardt\n",
    "    try:\n",
    "        print(f\"\\n{'='*20} LEVENBERG-MARQUARDT {'='*20}\")\n",
    "        results['levenberg_marquardt'] = levenberg_marquardt_solve(V, L_rho)\n",
    "    except Exception as e:\n",
    "        print(f\"Levenberg-Marquardt failed: {e}\")\n",
    "        results['levenberg_marquardt'] = None\n",
    "\n",
    "    return results"
   ],
   "id": "ced1aa445092295b",
   "outputs": [],
   "execution_count": 7
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T03:12:16.973908Z",
     "start_time": "2025-07-31T03:12:13.279279Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# 使用示例\n",
    "    # 加载数据\n",
    "S = np.load(path+'S_matrix.npy').astype(np.complex128)\n",
    "b = np.load(path+'b.npy').astype(np.complex128)\n",
    "V = np.load(path+'V.npy').astype(np.complex128)\n",
    "L_rho = np.load(path+'L_rho.npy').astype(np.complex128)\n",
    "\n",
    "    # 运行所有方法比较\n",
    "results = compare_all_methods(V, L_rho)\n",
    "\n",
    "    # 或者单独使用最推荐的方法\n",
    "best_solution = high_precision_solve(V, L_rho, method='svd_truncated', tolerance=1e-15)\n",
    "pass"
   ],
   "id": "fdfd2c9bbf4da7da",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "============================================================\n",
      "COMPARING ALL HIGH-PRECISION METHODS\n",
      "============================================================\n",
      "\n",
      "==================== SVD_TRUNCATED ====================\n",
      "rank(S) = 2\n",
      "rank([S | b]) = 2\n",
      "Condition number of S: 9.7122e+80\n",
      "\n",
      "=== SVD Truncated Solve ===\n",
      "Total singular values: 273\n",
      "Valid singular values (> 1.00e-12 * σ_max): 1\n",
      "Smallest valid singular value: 2.0510e+09\n",
      "Largest singular value: 2.0510e+09\n",
      "Condition number (before truncation): 1.4153e+16\n",
      "Residual norm for S @ a - b: 3.4849e-02\n",
      "\n",
      "==================== QR_PIVOTING ====================\n",
      "rank(S) = 2\n",
      "rank([S | b]) = 2\n",
      "Condition number of S: 9.7122e+80\n",
      "\n",
      "=== QR Pivoting Solve ===\n",
      "Matrix rank from QR: 2\n",
      "Smallest valid diagonal element: 2.3202e-02\n",
      "Largest diagonal element: 3.7256e+04\n",
      "Method qr_pivoting failed: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 4 is different from 273)\n",
      "\n",
      "==================== ITERATIVE_REFINEMENT ====================\n",
      "rank(S) = 2\n",
      "rank([S | b]) = 2\n",
      "Condition number of S: 9.7122e+80\n",
      "\n",
      "=== Iterative Refinement Solve ===\n",
      "Initial solution computed\n",
      "Iteration 1: Residual norm = 1.4462e-10\n",
      "Iteration 2: Residual norm = 2.0580e-11\n",
      "Iteration 3: Residual norm = 8.2318e-11\n",
      "Iteration 4: Residual norm = 2.0580e-11\n",
      "Iteration 5: Residual norm = 2.0580e-11\n",
      "Iteration 6: Residual norm = 2.0580e-11\n",
      "Iteration 7: Residual norm = 2.1084e-16\n",
      "Final residual norm for S @ a - b: 2.1084e-16\n",
      "\n",
      "==================== TIKHONOV_ADAPTIVE ====================\n",
      "rank(S) = 2\n",
      "rank([S | b]) = 2\n",
      "Condition number of S: 9.7122e+80\n",
      "\n",
      "=== Adaptive Tikhonov Regularization ===\n",
      "Optimal regularization parameter: 7.1969e-07\n",
      "Residual norm for S @ a - b: 1.0603e-06\n",
      "Warning: Solution has significant imaginary components\n",
      "Max imaginary part: 2.4075e-09\n",
      "\n",
      "==================== MULTIPLE_PRECISION ====================\n",
      "rank(S) = 2\n",
      "rank([S | b]) = 2\n",
      "Condition number of S: 9.7122e+80\n",
      "\n",
      "=== Multiple Precision Solve ===\n",
      "Using precision: complex128\n",
      "Residual norm for S @ a - b: 1.0604e-06\n",
      "\n",
      "==================== LEVENBERG-MARQUARDT ====================\n",
      "\n",
      "=== Levenberg-Marquardt Solve ===\n",
      "Levenberg-Marquardt failed: name 'ComplexWarning' is not defined\n",
      "rank(S) = 2\n",
      "rank([S | b]) = 2\n",
      "Condition number of S: 9.7122e+80\n",
      "\n",
      "=== SVD Truncated Solve ===\n",
      "Total singular values: 273\n",
      "Valid singular values (> 1.00e-15 * σ_max): 2\n",
      "Smallest valid singular value: 9.7951e-04\n",
      "Largest singular value: 2.0510e+09\n",
      "Condition number (before truncation): 1.4153e+16\n",
      "Condition number (after truncation): 2.0939e+12\n",
      "Residual norm for S @ a - b: 1.0604e-06\n"
     ]
    }
   ],
   "execution_count": 8
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}
